home *** CD-ROM | disk | FTP | other *** search
/ Trading on the Edge / Trading On The Edge - CD-ROM Toolkit (Wayzata Technology)(2031)(1994).bin / pc / shared / economet.m < prev    next >
Text File  |  1994-01-05  |  16KB  |  617 lines
  1. (*:Version: Mathematica 1.2 *)
  2.  
  3. (*:Name: Econometrics *)
  4.  
  5. (*:ID:            *)
  6.  
  7. (*:Title: Econometrics.m *)
  8.  
  9. (*:Author: David A. Belsley, Department of Economics, Boston College and Center for Computational Research in Economics and Management Science, MIT, 1990.*)
  10.  
  11. (*:Legal:          *)
  12.  
  13. (*:Reference:      *)
  14.  
  15. (* Copyright, David A. Belsley, 1989. Version 1.0. *)
  16.  
  17. BeginPackage["Econometrics`"]
  18. Reg::usage = "Reg[y_,x_,opts___] runs a regression of y on
  19. x.  y is a vector, x is a matrix or a list of vectors,
  20. matrices or lists of vectors.  Include \"const\" (quotes
  21. are essential here) in x list if there is to be an
  22. intercept.  y and x may contain missing observations
  23. (any non-numeric entries).\r \r
  24. Options[Reg] = {displayOutput -> True, varDCom ->False,
  25. returnValues -> Level1}"
  26.  
  27. TwoStage::usage = "TwoStage[y_,yincluded_,xincluded_,
  28. xexlcuded_,opts___] does two stage least squares.  y must
  29. be a vector (or 1xn matrix)."
  30.  
  31. IV::usage = "IV[y_,includedVariables_,instrumentalVariables_,
  32. opts___] does IV estimation.  If number of instVariables
  33. equals yincluded, instVariables become instruments, else
  34. the instruments are the yhats of includedy regressed on
  35. includedx's and instVariables. y must be vector or 1xn
  36. matrix."
  37.  
  38. lag::usage = "lag[list_,n_,m_] creates lags (leads) of list.
  39. Negative values denote leads. lag[list] is list lagged once.
  40. lag[list,n] is list lagged n times. lag[list,n,m] creates a
  41. matrix of list lagged n through m>n times.  This version
  42. preserves the length of the data series, adding 'na's as needed."
  43.  
  44. Mean::usage = "Mean[list_List], the arithmetic mean of list."
  45.  
  46. Correlation::usage = "Correlation[x1_List, x2_List], gives the
  47. correlation between any two lists of the same length."
  48.  
  49.  
  50.  
  51. (* Aliases make simpler names globally accessable  *)
  52. displayOutput ::= Econometrics`private`displayOutput
  53. varDCom ::= Econometrics`private`varDCom
  54. returnValues ::= Econometrics`private`returnValues
  55.  
  56. Begin["`private`"]
  57.  
  58. Reg[y_,x_,opts___] :=
  59. Block[{beta,sig2,sigEE,xx,se,n,g,p,missObs,err,yhat,
  60.     ydata,xdata,yCon,xCon,retainedCases,temp},
  61.     
  62.     If[TrueQ[!(temp=doRegData[y,x])],
  63.         Return[False],
  64.         {ydata,xdata,n,p,g,yCon,xCon,missObs,retainedCases} =
  65.                 temp
  66.     ];
  67.     
  68.     {beta,err,yhat,sig2,se,xx} = doRegCalculations[ydata,xdata,
  69.                             n,p,g,missObs,retainedCases,opts];
  70.  
  71.     If[TrueQ[missObs!=0],
  72.         err = Map[restore[#,retainedCases]&,err];
  73.         yhat = Map[restore[#,retainedCases]&,yhat] 
  74.     ];
  75.                                 
  76.     doReturn[opts]    (* must not have semicolon *)
  77. ]
  78.  
  79. Reg::vecx = "independent variables are not of proper form"
  80. Reg::vecy = "dependent variable(s) not of proper form"
  81. Reg::doff = "inadequate degrees of freedom"
  82. Reg::con2 = "Constant used twice";
  83.  
  84. Options[Reg] := {displayOutput -> True, varDCom -> False,
  85.     returnValues -> Level1};
  86.     
  87. Options[SimEqns] := {estimator -> TwoStage}
  88. doReturn[opts___] :=
  89.     Switch[ToString[returnValues/.{opts}/.Options[Reg]],
  90.         "Level2",
  91.             Return[Map[reduceRank,{beta,err,yhat,sig2,se,xx}]],
  92.         _,
  93.             Return[Map[reduceRank,{beta,err,yhat,sig2}]]
  94.     ]
  95.  
  96. doRegCalculations[yMat_,xMat_,n_,p_,g_,missObs_,    
  97.                             retainedCases_,opts___] :=
  98. Block[{v,sv,yhat,beta,xx,
  99.         err,sigEE,sig2,se,rsqr,rbar2,runcen},
  100.  
  101.     {v,sv,beta,xx} = doReg[yMat,xMat,opts];
  102.     
  103.     {yhat,err,sigEE,sig2,se,rsqr,rbar2,runcen} = 
  104.         doStats[yMat,xMat,beta,xx,n,p,g];
  105.     
  106.     If[displayOutput/.{opts}/.Options[Reg],
  107.         displayOut[beta,se,err,sig2,xCon,n,p,g,missObs,
  108.             retainedCases,rsqr,rbar2,runcen]
  109.     ];  
  110.     
  111.     If[varDCom/.{opts}/.Options[Reg],
  112.         doVarDCom[v,sv,p]
  113.     ];
  114.     
  115.     Return[{beta,err,yhat,sig2,se,xx}]
  116.         (* Note: all values returned here are matrices -- even
  117.             if row matrices -- not vectors. reduceRank is used
  118.             at doReturn to reduce ranks if needed. *)
  119. ]
  120.  
  121. doReg[yMat_,xMat_,opts___] :=
  122. Block[{u,sv,v,beta,xx,normx,S,
  123.             scaling = varDCom/.{opts}/.Options[Reg]},
  124.     
  125.     If[TrueQ[scaling],
  126.         normx = Map[EuclideanNorm,xMat];
  127.         S = DiagonalMatrix[1.0/normx]
  128.     ];
  129.     
  130.     {v,sv,u} = SingularValues[
  131.                     If[TrueQ[scaling],
  132.                         xMat/normx,
  133.                         xMat
  134.                     ],
  135.                 Tolerance -> 10.^-30];
  136.                 
  137.     u = Transpose[v].DiagonalMatrix[1/sv].u;(* pseudo inv *)
  138.     
  139.     beta = If[TrueQ[scaling],
  140.                 yMat.Transpose[u].S,
  141.                 yMat.Transpose[u]
  142.             ];
  143.             
  144.     xx = If[TrueQ[scaling],
  145.             S.Transpose[v].DiagonalMatrix[1/sv^2].v.S,
  146.             Transpose[v].DiagonalMatrix[1/sv^2].v
  147.         ];
  148.     Return[{v,sv,beta,xx}]
  149. ]
  150.  
  151. doStats[yMat_,xMat_,beta_,xx_,n_,p_,g_] :=
  152. Block[{i,temp,yhat,err,sigEE,sig2,se,rsqr,rbar2,runcen},
  153.     yhat = beta.xMat;
  154.     err = yMat-yhat;
  155.     sigEE = err.Transpose[err]/(n-p);        
  156.     sig2 = Table[sigEE[[i,i]],{i,g}];
  157.     temp = Sqrt[Table[xx[[i,i]],{i,p}]];
  158.     se = DiagonalMatrix[Sqrt[sig2]].Table[temp,{g}];
  159.     rsqr = Table[Correlation[yhat[[i]],yMat[[i]]]^2,{i,g}];
  160.     rbar2 = 1.0-((n-1)/(n-p))(1.0-rsqr);
  161.     runcen = Table[yhat[[i]].yhat[[i]]/yMat[[i]].yMat[[i]],{i,g}];
  162.  
  163.     Return[{yhat,err,sigEE,sig2,se,rsqr,rbar2,runcen}]
  164. ]
  165.  
  166. displayOut[beta_,se_,err_,sig2_,xCon_,n_,p_,g_,missObs_,
  167.         retainedCases_,rsqr_,rbar2_,runcen_] :=
  168. Block[{i,j,f,labels=Array[f,p],er = err},
  169.     For[i=1,i<=p,i++,
  170.         labels[[i]] = StringJoin["X[",ToString[i],"]"]
  171.     ];
  172.     If[xCon,labels[[1]] = "Const"];
  173.     For[i=1,i<=g,i++,
  174.         Print[StringForm["Dependent variable ``",i] ];
  175.         If[TrueQ[missObs!=0],
  176.             er = Map[restore[#,retainedCases]&,er];
  177.             dw = Apply[Plus,Select[(Drop[er[[i]],1]-
  178.                 Drop[er[[i]],-1])^2,NumberQ]]/Apply[Plus,
  179.                 Select[er[[i]]^2,NumberQ]],
  180.             dw = Apply[Plus,(Drop[er[[i]],1]-
  181.                 Drop[er[[i]],-1])^2]/((n-p)*sig2[[i]])
  182.         ];
  183.     
  184.         Print[StringForm["RSquared = ``  RBarSquared = ``",
  185.             rsqr[[i]],rbar2[[i]]]];
  186.         Print[StringForm["R2uncentered = ``  SER = ``",
  187.             runcen[[i]],Sqrt[sig2[[i]]] ]];
  188.         Print[StringForm["Num of Observations = ``
  189.                 Degrees of Freedom = ``",n,n-p] ];
  190.         Print[StringForm["dw = `` with `` missing obs.",
  191.             dw, missObs] ];
  192.         Print[" "];
  193.         Print[StringForm["var \t\tcoef.\t\tst. err.\t\tt"]];
  194.         Print[" "];
  195.         
  196.         For[j=1,j<=p,j++,
  197.             Print[StringForm["``\t\t``\t\t``\t\t``",
  198.                 labels[[j]], N[beta[[i,j]],5],
  199.                     N[se[[i,j]],5],
  200.                     N[beta[[i,j]]/se[[i,j]],5] ] ] ];
  201.         Print[" "];
  202.         Print[" "]
  203.     ]
  204. ]
  205.  
  206. doVarDCom[v_,sv_,p_] :=
  207. Block[{vMat,i,j,temp,tempStr},    
  208.  
  209.     vMat = (Transpose[v].DiagonalMatrix[1.0/sv])^2;
  210.     vMat = Map[(#/Apply[Plus,#])&,vMat]; 
  211.     
  212.     
  213.     
  214. (* Form condition indexes and output matrix with
  215.     condition indexes in first row and the variance-
  216.     decomposition proportions in the next p rows *)
  217.                  
  218.      vMat = Truncate[vMat,3];
  219.          PrependTo[vMat,Max[sv]/sv];
  220.          vMat = Transpose[Reverse[Sort[Transpose[vMat]]]];
  221.  
  222.     Print[" "];
  223.     Print["Variance-decomposition proportions"];
  224.     Print[" "];
  225.     Print["                  Condition indexes"];
  226.     Print[" "];
  227.     
  228.     
  229.     tempStr = "    ";
  230.     For[i=1,i<=p,i++,
  231.         tempStr = StringJoin[tempStr,"\t",ToString[FortranForm
  232.             [N[vMat[[1,i]],3]]]]
  233.     ];
  234.     Print[tempStr];
  235.     
  236.     Print["  "];
  237.     For[i=1,i<=p,i++,
  238.         tempStr = ToString[StringForm["V[``]",i]];
  239.         For[j=1,j<=p,j++,
  240.             tempStr = StringJoin[tempStr,"\t",ToString[StringForm["``",
  241.                 If[(temp = vMat[[i+1,j]])<1.0 && temp > 0.0, temp, 
  242.                 If[temp == 0.0, "0.000","1.000"]]] ]]
  243.         ];
  244.         Print[tempStr];
  245.     ];
  246. ]
  247.  
  248. doRegData[yblock_,xblock_] :=
  249. Block[{ydata,xdata,n,p,g,xCon,yCon,missObs,retainedCases,
  250.     temp,ny,nx},
  251.     If[TrueQ[!(temp = checkInput[yblock])],
  252.         Message[MessageName[Reg,"vecy"]];
  253.         Return[False],
  254.         {ydata,ny,g,yCon} = temp
  255.     ];
  256.     
  257.     If[TrueQ[!(temp = checkInput[xblock])],
  258.         Message[MessageName[Reg,"vecx"]];
  259.         Return[False],
  260.         {xdata,nx,p,xCon} = temp
  261.     ];
  262.     
  263.     If[yCon && xCon,
  264.         Message[MessageName[Reg,"con2"]];
  265.         Return[False]
  266.     ];
  267.  
  268.     n = Min[nx,ny];
  269.     
  270.     If[yCon,PrependTo[ydata,Table[1.0, {n}]] ];
  271.     If[xCon,PrependTo[xdata,Table[1.0, {n}]] ];
  272.  
  273.     If[TrueQ[!(temp=missingObservations[Join[ydata,xdata],n,p])],
  274.         Message[MessageName[Reg,"doff"]];
  275.         Return[False],
  276.         ydata = Take[temp[[1]],{1,g}];
  277.         xdata = Drop[temp[[1]],{1,g}];
  278.         {retainedCases,missObs,n} = Drop[temp,1];
  279.     ];
  280.     
  281.     Return[{ydata,xdata,n,p,g,yCon,xCon,missObs,retainedCases}]
  282. ]
  283.  
  284. checkInput[dataBlock_] :=
  285.  
  286. (*  Produces a data matrix from an input of any combination 
  287. of lists, matrices or irregular arrays.  Sets globals for 
  288. checkParts[] which is called recursively.  On return, data is 
  289. the possibly irregular data array, nobs is the length of the
  290. shortest row in data, varCount is the number of rows, and Con 
  291. is set True if a "const" is encountered and otherwise False *)
  292.  
  293. Block[{varCount=0,data={},Con=False,nobs=Infinity},
  294.     If[!checkParts[dataBlock], Return[False] ];
  295.             
  296.     Return[{data,nobs,varCount,Con}]
  297. ]
  298.  
  299. checkParts[el_] :=
  300.  
  301. (*  Recursive routine to coalesce divergent input forms into
  302. a single, possibly irregular array.  Globals set by calling
  303. routine: Con,data,nobs, and varCount. nobs is set to minimum 
  304. of all data lengths.  Con is set True if a "const" is included 
  305. among el. Returns True unless an empty series is encountered.*)
  306.  
  307. Block[{nt,pt,i},
  308.  
  309.     If[MemberQ[{"const","Const","constant","Constant"},el],
  310.         If[Con,
  311.             Message[MessageName[Reg,"con2"]];
  312.             Return[False],
  313.             varCount++;
  314.             Con = True;
  315.             Return[True]
  316.         ]
  317.     ];
  318.     If[TrueQ[Depth[el] == 1],
  319.         Return[False]  
  320.     ];
  321.     If[VectorQ[el],
  322.         nobs = Min[nobs,Length[el]];
  323.         varCount++;
  324.         AppendTo[data,el];
  325.         Return[True]
  326.     ];
  327.     If[MatrixQ[el]&&(Depth[el]==3), 
  328. (* This second condition is necessary only because
  329. MatrixQ[] works incorrectly  *)
  330.         {nt,pt} = Dimensions[el];
  331.         If[nt < pt,
  332.             {nt,pt} = {pt,nt};
  333.             data = Join[data,el],
  334.             data = Join[data,Transpose[el]]
  335.         ];
  336.         nobs = Min[nobs, nt];
  337.         varCount += pt;
  338.         Return[True]
  339.     ]; 
  340.     For[i=1,i<=Length[el],i++,
  341.         If[!checkParts[el[[i]]],
  342.             Return[False]
  343.         ]
  344.     ];
  345.     Return[True]
  346. ]
  347.  
  348. missingObservations[dataBlock_,blockLength_,testLength_] :=
  349.  
  350. (* routine to remove incomplete observations - ones with
  351. non-numerical entries.  Globals: n, p, xdata, ydata, 
  352. retainedCases,missObs. At end, xdata and ydata contain data 
  353. suitable for regression, and caseRecord contains indexes of
  354. retained observations. Returns True unless degrees of freedom 
  355. become zero or below. *)
  356.  
  357. Block[{i=1,j,tempData,missingCount=0,nobs=blockLength,
  358.                 retainedCases},
  359.     tempData = Map[Take[#,nobs]&,dataBlock]//N;
  360.     PrependTo[tempData,Table[j,{j,nobs}]];
  361.     tempData = Transpose[tempData];
  362.     
  363.     While[i<=nobs,
  364.         If[TrueQ[Scan[If[!NumberQ[#], 
  365.                 Return[False]]&,tempData[[i]] ]==False],
  366.             tempData = Drop[tempData,{i,i}];
  367.             missingCount++;
  368.             nobs--,
  369.             i++
  370.         ];
  371.         If[nobs<=testLength,
  372.             Return[False]
  373.         ]            
  374.     ];
  375.     tempData = Transpose[tempData];
  376.         
  377.     retainedCases = tempData[[1]];
  378.     tempData = Drop[tempData,1];
  379.     Return[{tempData,retainedCases,missingCount,nobs}]
  380. ] 
  381.  
  382. restore[series_,index_] :=
  383. Block[{temp,i},
  384.     temp = Table["na",{index[[-1]]}];
  385.     For[i=1,i<=Length[index],i++,
  386.         temp[[index[[i]]]] = series[[i]]
  387.     ];
  388.     Return[temp]
  389. ]
  390.  
  391. IV[y_,includedVariates_,instrumentalVariables_,opts___]:=
  392.     TwoStage[y,includedVariates,{},instrumentalVariables,opts,
  393.         estimator -> IV]
  394. TwoStage[y_,includedy_,includedx_,excludedx_,opts___] :=
  395. Block[{yhat,beta,xx,ydata,y1data,x1data,x2data,n,g,g1,k1,k2,
  396.     Con,retainedCases,err,missObs,sig2,se,sigEE,rsqr,rbar2,
  397.     runcen,temp},
  398.     
  399.     (* Prepare data *)
  400.     If[TrueQ[!(temp = do2SData[y,includedy,includedx,
  401.                                             excludedx,opts])],
  402.         Return[False],
  403.         {ydata,y1data,x1data,x2data,n,g,g1,k1,k2,Con,
  404.             retainedCases,missObs} = temp;
  405.     ];
  406.     
  407.     Switch[estimator/.{opts}/.Options[SimEqns],
  408.         IV,
  409.             If[k2>g1,
  410.                 temp = doReg[y1data,x2data];
  411.                 x2data = temp[[3]].x2data
  412.             ];
  413.             {beta,xx} = doIV[ydata,y1data,x2data],
  414.             
  415.         _,  (* Default is TwoStage  *)
  416.                 (* form first stage *)
  417.             temp = doReg[y1data,Join[x1data,x2data]];
  418.             yhat = temp[[3]].Join[x1data,x2data];
  419.             
  420.                 (* Adjust constant to first row if present *)
  421.             If[Con,
  422.                 PrependTo[yhat,x1data[[1]]];
  423.                 PrependTo[y1data,x1data[[1]]];
  424.                 x1data = Drop[x1data,1];
  425.             ];
  426.                 (* Form estimator -- or second stage  *)
  427.             temp = doReg[ydata,Join[yhat,x1data]];
  428.             beta = temp[[3]];
  429.             xx = temp[[4]]
  430.     ];
  431.     
  432.  
  433.     (* second-stage statistics *)
  434.     {yhat,err,sigEE,sig2,se,rsqr,rbar2,runcen} = 
  435.         doStats[ydata,Join[y1data,x1data],beta,xx,n,g1+k1,g];
  436.         
  437.     If[displayOutput/.{opts}/.Options[Reg],
  438.         displayOut[beta,se,err,sig2,Con,n,g1+k1,g,missObs,
  439.             retainedCases,rsqr,rbar2,runcen]
  440.     ];
  441.  
  442.     If[varDCom/.{opts}/.Options[Reg],
  443.         doVarDComIV[xx,reduceRank[se],reduceRank[sig2],g1+k1]
  444.     ];
  445.     
  446.     If[TrueQ[missObs!=0],
  447.         err = Map[restore[#,retainedCases]&,err];
  448.         yhat = Map[restore[#,retainedCases]&,yhat] 
  449.     ];
  450.     
  451.     doReturn[opts] (* must not have semicolon *)
  452. ]
  453.  
  454. doVarDComIV[xx_,se_,sig2_,p_] :=
  455. Block[{u,sv,v,S,vMat},
  456.     S = DiagonalMatrix[1.0/se];
  457.     vMat = S.(sig2*xx).S;
  458.     {v,sv,u} = SingularValues[vMat,Tolerance -> 10.^-30];
  459.     doVarDCom[v,Sqrt[sv],p]
  460. ]
  461.  
  462. do2SData[y_,y1_,x1_,x2_,opts___] :=
  463. Block[{ydata,y1data,x1data,x2data,n,g,g1,k1,k2,Con,
  464.     x1Con,x2Con,y1Con,temp,ny1,nx1,nx2,retainedCases,missObs},
  465.     
  466.     If[TrueQ[!(temp = checkInput[y])],
  467.         Message[MessageName[TwoStage,"ydat"]];
  468.         Return[False],
  469.         {ydata,n,g,Con} = temp;
  470.     ];
  471.     
  472.     If[g!=1,
  473.         Message[MessageName[TwoStage,"ydat"]];
  474.         Return[False]
  475.     ];
  476.     
  477.     If[TrueQ[!(temp = checkInput[y1])],
  478.         Message[MessageName[TwoStage,"y1dat"]];
  479.         Return[False],
  480.         {y1data,ny1,g1,y1Con} = temp;
  481.     ];
  482.     
  483.     If[TrueQ[!(temp = checkInput[x2])],
  484.         Message[MessageName[TwoStage,"x2dat"]];
  485.         Return[False],
  486.         {x2data,nx2,k2,x2Con} = temp;
  487.     ];
  488.     
  489.     If[k2<g1,
  490.         Message[MessageName[TwoStage,"unid"]];
  491.         Return[False]
  492.     ];
  493.     Switch[estimator/.{opts}/.Options[SimEqns],
  494.         IV,
  495.             x1data = {};
  496.             nx1 = Infinity;
  497.             k1 = 0;
  498.             x1Con = False;
  499.             Con = y1Con,
  500.         _,
  501.             If[TrueQ[!(temp = checkInput[x1])],
  502.                 Message[MessageName[TwoStage,"x1dat"]];
  503.                 Return[False],
  504.                 {x1data,nx1,k1,x1Con} = temp;
  505.             ];
  506.             If[(x1Con&&x2Con)||(y1Con&&x1Con)||(y1Con&&x2Con),
  507.                 Message[MessageName[TwoStage,"con2"]];
  508.                 Return[False]
  509.             ];
  510.             Con = x1Con
  511.     ];
  512.  
  513.     n = Min[n,ny1,nx1,nx2];
  514.     
  515.     If[y1Con,PrependTo[y1data,Table[1.0, {n}]] ];
  516.     If[x1Con,PrependTo[x1data,Table[1.0, {n}]] ];
  517.     If[x2Con,PrependTo[x2data,Table[1.0, {n}]] ];
  518.     
  519.     If[TrueQ[!(temp = missingObservations[Join[ydata,
  520.         y1data,x1data,x2data],n,k1+k2])],
  521.             Message[MessageName[TwoStage,"doff"]];
  522.             Return[False],
  523.             {temp,retainedCases,missObs,n} = temp;
  524.     ];
  525.     
  526.     ydata = Take[temp, {1,1}];
  527.     y1data = Take[temp, {2,g1+1}];
  528.     x1data = Take[temp, {g1+2,g1+1+k1}];
  529.     x2data = Take[temp, {-k2,-1}];
  530.     
  531.     Return[{ydata,y1data,x1data,x2data,n,g,g1,k1,k2,Con,
  532.         retainedCases, missObs}]
  533. ]
  534.  
  535. doIV[y_,z_,w_] :=
  536. Block[{delta,xx,zwInv},
  537.     zwInv = Inverse[z.Transpose[w]];
  538.     xx = Transpose[zwInv].
  539.             (w.Transpose[w]).zwInv;
  540.     delta = (y.Transpose[w]).zwInv;
  541.     Return[{delta,xx}]
  542. ]
  543.  
  544. TwoStage::ydat := "y variate incorrect."
  545. TwoStage::y1dat := "Included variables or included
  546. endogenous variates incorrect."
  547. TwoStage::x1dat := "Included exogenous variates incorrect."
  548. TwoStage::x2dat := "Instrumental variables or excluded
  549. exogenous variates incorrect."
  550. TwoStage::unid := "Inadequate instruments -- underidentified."
  551. TwoStage::con2 := "Constant used twice."
  552. TwoStage::doff := "Inadequate degrees of freedom."
  553.  
  554. (* Utility Programs for Reg *)
  555.  
  556. Mean[list_List] := N[Apply[Plus,list]/Length[list]]
  557.  
  558. Correlation[x1_List,x2_List] :=
  559. Block[{x1bar,x2bar},
  560.     x1bar = Mean[x1];
  561.     x2bar = Mean[x2];
  562.     (x1-x1bar).(x2-x2bar)/Sqrt[(x1-x1bar).(x1-x1bar)*
  563.                                     (x2-x2bar).(x2-x2bar)]
  564. ]/; Length[x1] == Length[x2]
  565.  
  566. EuclideanNorm[list_] := N[Sqrt[list.list]]
  567.  
  568. printString[str_,args__] := WriteString["stdout",StringForm[
  569.     str, args] ]
  570.  
  571. Truncate[x_,n_] := N[Round[x 10.^n] 10.^-n]
  572.  
  573. reduceRank[x_] := If[MatrixQ[x]&&(Min[Dimensions[x]]==1)||
  574.     VectorQ[x]&&(Length[x]==1),
  575.     Return[x[[1]]],Return[x] ]
  576.     
  577. (* lag[list_List,n_:1,m_:undefined] :=
  578. Block[{i,temp},
  579.     If[TrueQ[m==undefined],
  580.         Which[
  581.             n<0,
  582.                 temp = Nest[Drop[#,1]&,list,-n],
  583.             n==0,
  584.                 temp = list,
  585.             n>0,
  586.                 temp = Drop[Nest[Prepend[#,"a"]&,list,n],-n]
  587.         ],
  588.         temp = Table[lag[list,i], {i,n,m}]
  589.     ];
  590.     temp = temp  (* This kludge corrects for Mathematica bug
  591.                 when Return[] is used in a function with
  592.                 optional parameters with defaults.  Remove
  593.                 temps and replace with Returns when bug
  594.                 is fixed. *)
  595. ]/;(m==undefined || m>n)  old version  *)
  596.  
  597. lag[list_List,n_:1,m_:undefined] :=
  598. Block[{i,temp},
  599.     If[TrueQ[m==undefined],
  600.         Which[
  601.             n<0,
  602.                 temp = Join[Drop[list,-n],Table["a",{-n}]],
  603.             n==0,
  604.                 temp = list,
  605.             n>0,
  606.                 temp = Join[Table["a",{n}],Drop[list,-n]]
  607.         ],
  608.         temp = Table[lag[list,i], {i,n,m}]
  609.     ];
  610.     Return[temp]
  611. ]/;(m==undefined || m>n)
  612.  
  613. End[]
  614.  
  615.  
  616. EndPackage[]
  617.